Collective stochastic resonance in shear-induced melting of sliding bilayers 
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The far-from-equilibrium dynamics of two crystalline two-dimensional monolayers driven past 
each other is studied using Brownian dynamics simulations. While at very high and low driving 
rates the layers slide past one another retaining their crystalline order, for intermediate range of 
drives the system alternates irregularly between the crystalline and fluid-like phases. A dynamical 
phase diagram in the space of interlayer coupling and drive is obtained. A qualitative understanding 
of this stochastic alternation between the liquid-like and crystalline phases is proposed in terms of 
a reduced model within which it can be understood as a stochastic resonance for the dynamics of 
collective order parameter variables. This remarkable example of stochastic resonance in a spatially 
extended system should be seen in experiments which we propose in the paper. 

PACS numbers: 82.70.Dd,05.45.Xt,62.20.Fe,81.40.Pq 



I. INTRODUCTION AND RESULTS 



A. Background 



in detail below, are quite distinct from the well-known 
stick-slip effect in atomically thin fluid films subjected to 
shear Ha HQ 111. 
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The shear flow of a solid is one of the most important 
and widely studied 0, 0, 0, 13 nonequilibrium phenom- 
ena in materials science, with relevance to such practi- 
cal problems as the yielding of materials, solid friction 
and even the mechanical properties of the earth's crust. 
Such flow takes place when solids are subjected to stresses 
which range from a few percent of the shear modulus to, 
in some cases, a value of the order of the shear modulus 
itself. It is particularly convenient to study such phe- 
nomena using very soft solids, where the desired stress 
to modulus ratio is easily achieved. Indeed, such studies 
open up new regimes in the physics of driven systems. 
A variety of such unconventional, ultra-soft solids have 
been studied, including packings of multilamellar vesi- 
cles 0, vortex lattices in type II superconductors @, 
and crystalline arrays, electrostatically or sterically sta- 
bilized, of colloidal particles in aqueous suspensions 0. 
Experiments on suspensions of interacting colloidal par- 
ticles under shear are of particular interest to us here, 
for the rich range of interesting phenomena they reveal, 
including the shear-induced distortion of the static struc- 
ture factor in the fluid state, and stick-sli p dy namics [8|, 
hysteresis 9] and shear induced melting |lft Hlj , in the 
crystalline state. It is likely that the properties of sheared 
crystals, as observed in macroscopic three-dimensional 
scattering studies or in time- or frequency-domain me- 
chanical measurements, are the average result of many in- 
termittent, spatially inhomogeneous internal events. Ac- 
cordingly, this paper focuses on such events, at the level 
of the relative motion of an adjacent pair of layers, since 
we believe that knowledge of these events will greatly aid 
our understanding of the mechanisms underlying phe- 
nomena such as shear-melting. We emphasize at the 
outset, to avert any confusion on this score, that the phe- 
nomena which our study uncovers, and which we discuss 



One popular approach to the study of sheared solids 
has been to consider an ordered layer (the adsorbate), 
dragged over a fixed, rigid periodic potential (the sub- 
strate), the latter representing an adjacent layer [TrjlTR 
Hs| . This description is clearly limited in its applicability 
since it rules out deformation of the substrate, although 
it is a reasonable starting point for experimental situa- 
tions in which the overlayer is much softer than the sub- 
strate. It is natural, and more general, to ask instead 
what happens when both adsorbate and substrate are 
dynamical, and organize themselves into various struc- 
tures, depending on interaction strengths, temperature 
and driving force, and it is in this spirit that our model 
is formulated. The case where both layers are compara- 
bly deformable, in particular, is clearly of relevance to 
sheared crystals. In all cases, each layer confronts a pe- 
riodic potential produced by the other layer, but both 
amplitude and phase of this periodic potential are dy- 
namical and change as a result of interactions, noise and 
driving force, giving rise to some remarkable collective ef- 
fects, reported briefly earlier [lj| and discussed in detail 
in this paper. 



Although the primary motivation for this paper was 
the problem of sheared colloidal crystals, there are two 
other classes of problem to which our study has a natural 
connection. One is the phenomenon of lane formation 
in counter-driven interacting particles po|. the other is 
the equilibrium modulated-liquid to solid transition of 
interacting particles in an external periodic potential. We 
will touch upon the relation of these problems to our work 
later in this paper. 
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B. Summary of models and results 

We report two detailed studies in this paper: first, a 
Brownian dynamics simulations of a many-particle model 
[T^| , henceforth referred to as the particle model, and sec- 
ond, a reduced model, introduced to get insight into the 
results of the particle model, consisting of just two de- 
grees of freedom [2l] , an order parameter amplitude and 
a strain field. The particle model consists of two species 
of particles, A and B, driven by a force F, of constant 
magnitude, in opposite directions say, along +x and —x 
respectively (Fig. . The Vaa and Vb b interactions are 
identical. The Vab interaction has the same form but 
is smaller by a factor e. This factor, in a phenomeno- 
logical way, incorporates the physics of the third direc- 
tion (see below) . All pairwise interactions are of screened 
Coulomb form, with the screening parameter so chosen 
that, when F = 0, each species in the absence of the other 
settles down in a macroscopically ordered triangular lat- 
tice configuration. The dynamics of the system is mod- 
eled by the overdamped Langevin equation for relatively 
sheared sets of particles and is monitored for different 
F and e. With the inter-layer coupling strength e held 
constant, on increasing the drive, we observe an interest- 
ing sequence of non-equilibrium states, namely, a sliding 
crystalline ordered state (Fig. |2J), a sliding melt-freeze 
state (characterized by alternate states of order and dis- 
order in time), followed again by a sliding ordered state. 
In the intermediate "melt-freeze" regime, for fixed drive, 
the residence time of the system in the ordered state de- 
creases and that of the disordered state increases as a 
function of e (Figs.|30J[SJ). The allowed nonequilibrium 
states are best understood in terms of a dynamic phase 
diagram of these states. We present such a dynamical 
nonequilibrium phase diagram (Fig. I10|l demarcating the 
three regimes (i) lower smooth sliding, (ii) alternating 
melt-freeze state, and (iii) upper smooth sliding state. 
The melt-freeze alternations are most pronounced in a 
window of driving force F and inter-layer coupling e val- 
ues. These melt-freeze cycles are strongly reminiscent of 
the time series of a system undergoing stochastic reso- 
nance |22l 1231 124. l25| and, to explore this aspect in more 
detail following ref . [2j| , we introduce the reduced model. 
Using the reduced model, we study the time evolution of 
the system using coupled time dependent Ginzburg Lan- 
dau equations for the order parameter and strain fields, 
as a function of a coupling parameter (a) entering the 
model equations and a drive (fi) analogous to e and F 
respectively. For a certain range of values of a, keeping a 
fixed, as a function of the drive parameter O, we observe 
three regimes (Fig. II II and Fig. ll2|) - a crystalline state ( 
nonzero order parameter value) , a bistable regime where 
the system alternates between the crystalline and liquid 
state (order parameter values being zero), followed again 
by a crystalline state. Keeping Q fixed at an optimum 
value, we find that the ratio of the average lifetime of the 
crystalline state to that of the the liquid state in the inter- 
mediate regime of bistability decreases as a is increased 



(Fig. 1131 and Fig. ll4|l . These observations are remarkably 
similar to the phenomenon observed in the particle model 
and indeed the phase diagrams of the two models fFig.HUI 
and Fig. I17f) correspond surprisingly well. Further, the 
reduced model exhibits a maximum in the signal to noise 
ratio at optimum values of the noise intensity (Fig. I16fl . 
thereby making the connection to stochastic resonance 
concrete pi I2I I2I I25T] . 

The paper is organized as follows. The Brownian dy- 
namics simulations of the particle model are described in 
Section II A and the results discussed in detail in section 

II B. This is followed by physical arguments in support 
of the behaviors observed. The reduced model [21| is 
introduced in Section III A and its results discussed in 

III B. Finally in Section IV we provide a discussion of 
our results, suggest how our observations may be verified 
experimentally and outline directions of future research. 



II. BROWNIAN DYNAMICS SIMULATIONS OF 
TWO ADJACENT MONOLAYERS 

A. Particle Model 

We consider two sets A and B of Brownian particles in 
two spatial dimensions, driven in the +x and -x directions 
respectively by a constant driving force with magnitude 
F as shown in Fig. ^ Pairwise interactions between par- 
ticles are described by potentials VaaM, UBs(r) and 
Vab(i"). We choose a rectangular box of dimensions 
L = (V3/2) x 20£ and W = 201, where I = {2V3n )- 1/2 , 
tiq being the mean number density of either species. All 
quantities we use are in nondimensional form. Lengths 
are non-dimensionalized by I and time by r = £ 2 /D,D 
being the Brownian diffusivity. Energy is scaled by ksT 
and force by ksT/£, where T is the temperature and ks 
the Boltzmann constant. If the two close packed planes, 
with species A and B in distinct planes, are slid past each 
other, it is highly plausible that Vab will be weaker than 
Vaa and Vb b (since the effect of the interlayer interac- 
tion on the two-dimensional in-plane motion of the parti- 
cles, which is what Vab encodes, is substantially smaller 
than that due to the intralayer interaction) and that the 
strength of Vab relative to Vaa and Vbb can be varied 
by increasing or decreasing the normal confining pres- 
sure. We therefore choose, for this work dimensionless 
pair potentials of the screened Coulomb form 

V AA (r) = V BB {r) = z-*Vab{t) = {U /r) exp(-Kr), 

where Uq (= 1.75 x 10 4 ) and the screening parameter 
(k£ = 0.5 ) are so chosen that each species in the absence 
of the other and any external driving force settles down 
in a triangular lattice configuration (Fig. QJ. The par- 
ticle positions then evolve according to the overdamped 
Langevin equations 

V(t + St) = V(<) + St[F^ + f/(R M (i)) + V(i)],(2) 
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FIG. 1: Schematic diagram of the model 



where R^ 1 = {xi^,yi^) is the coordinate of the ith particle 
of the species /i (=4or5), and ±Fx (+ for A, - for B) 
is the external driving force on the ith particle of type \x. 

V(Ry = -^V^ !/3 -(R i -R^) (3) 

are the interparticle forces, and are Gaussian white 
noise sources with zero mean obeying the fluctuation dis- 
sipation relation which in nondimensional form reads 

{h fl i (0)hJ(t))=2\S flu S i ^S(t). (4) 

It is trivial to generalize to the case where species A and 
B differ but we have chosen them to be same in this 
model. We study the time evolution of this system as a 
function of the drive keeping e constant for several values 
of e. The dimcnsionless time step used in our simulations 
is St = 6.4 x 1(T 6 . 



B. Simulation Results 

The results reported in our study are generally for 
10 6 — 10 7 time steps after the initial transients ~ 10 4 
steps are discarded. Over this time, the A and B lat- 
tices sweep through each other a few to several hundred 
times depending upon the magnitude of the drive. In 
order to drift under the action of the driving force F, 
the particles have to overcome a barrier of the order of 
Vab{@) arising from interaction with the nearest neigh- 
bors of the opposite species. Thus, although F is itself 
dimensionless, it is appropriate to state the results in 
terms of the physically relevant dimensionless combina- 
tion Fd = F£/Vab(£)- However, for the phase diagram in 
the e-F variables, we have used the dimensionless combi- 
nation Fd* = F£€/Vab(£), as Vab already incorporates 
a factor of e in its definition. The structure and dynam- 
ics of the system has been monitored through snapshots 
of configurations, drift velocities Vd, particle-averaged lo- 
cal velocity variances ((Sv) 2 }, pair correlation functions 
gfiui^) as functions of separation r, and time dependent, 
but equal time structure factors S^^k, t) as functions of 
wavevector k (fj, and v range over A, B). Each point in 
the <7 M i/(r), (r = x, y) is an average over 100 data points, 




FIG. 2: Simulation images of macroscopically ordered lattices 
drifting through each other for e — 0.05, = 0.0438. 



recorded at times separated by 5Q6t. In the absence of 
the driving force ( i.e., at Fd = 0), the system is an 
imperfectly ordered crystal. The application of a small 
nonzero Fd, well below the apparent threshold for per- 
ceptible macroscopic relative motion of the two lattices 
(of A and B particles respectively), facilitates particles 
that are initially in unfavorable positions to re-organize 
and move to favorable locations leading to a small move- 
ment in these regions. After these transient motions the 
system settles down into a macroscopically ordered struc- 
ture with both components showing perfect long-range 
crystalline order sustained over distances of the order of 
the system size. There is no further relative drift of A 
and B except perhaps a tiny activated creep which we 
cannot resolve. Thereafter, keeping interaction strengths 
and temperature fixed, the driving force F displays three 
threshold values Fi, i = 1, 2, 3 corresponding to the lower 
bounds of three states -a lower sliding crystalline state, 
a melt-freeze state and an upper sliding crystalline state. 
The characteristic features of each of these three states 
are mentioned below. 

1. When Fd crosses the first threshold F\, the A and 
B components acquire a measurable, macroscopic rela- 
tive drift velocity Vd- The drift velocity shows a smooth 
change at this threshold value with the velocit y flu ctua- 
tions 2G] showing a pronounced enhancement 19] char- 
acteristic of depinning. (See Fig. 3 and 4 of Ref.[l9|.) 
This is likely to be a strong crossover rather than a true 
transition. Each particle faces a finite barrier to motion, 
so that at any nonzero temperature, particles can cross 
the barrier individually in an incoherent manner for ar- 
bitrarily small Fi. The barrier for creep velocity should 
thus be finite even in the limit of infinite system size. In 
the region F\ < F < F2, both A and B components are 
well-ordered, drifting crystals (Fig. 0). The two compo- 
nents slide smoothly past each other in lanes of width 
equal to the inter-particle distance, with negligible dis- 
tortion or disorder. We comment later on the connection 
to the lane formation work of [2Cj. 

2. In a window of driving forces F% < F < F3, we 
observe intriguing stochastic alternations of the system 
between an ordered and a disordered state which are the 
melt-freeze cycles (Figs. 0] El We have used snap- 
shots of configurations Figs. and |7| the pair correla- 
tion function g^ v {v) : Y = x, y (Figs. |H1 and |3J) and the 
equal-time but time-dependent structure factor S(k, t) 
to characterize these dynamical states of order and dis- 
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FIG. 3: The structure factor height (averaged over 1st ring of 
maxima) as a function of time in the melt-freeze cycle state 
for e = 0.05, F d * = 0.8167. 



FIG. 5: The structure factor height (averaged over 1st ring of 
maxima) as a function of time in the melt-freeze cycle state 
for e = 0.06, F d * = 1.4350. 
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FIG. 4: The structure factor height (averaged over 1st ring of 
maxima) as a function of time in the melt-freeze cycle state 
for e = 0.02, F d * = 0.8167. 



order. As the first two methods are representative of the 
instantaneous state of the system, for monitoring these 
cycles continuously we use the peak height of the (short- 
time averaged) static structure factor S(k,t). Indeed, 
the essential features of these cycles, namely the persis- 
tence duration, fluctuations in the extent of order etc, are 
best captured through the time dependence of 5(k, t) as 
we shall see. A typical such plot is shown in Fig. [3] for 
an optimum value of e = 0.05. It is clear that S(k,t) 
alternates between long stretches of crystal-like ( corre- 
sponding to large values of S(k, t)) and comparably long 
stretches of liquid-like values ( small values of 5(k, t) ) 
as the simulation progresses. This is strikingly different 
from stick-slip alternations, in which the melted (slip) 
state is considerably shorter than the crystalline (stick) 
state ^2 • We discuss this comparison later in section 
IV. For smaller values of e, even as the persistence of the 
crystalline state is enhanced, the extent of ordering itself 
is less pronounced as indicated by the lower values of S 
( compared to that for the crystalline phase correspond- 
ing to e = 0.05). Concomitantly the liquid- like state also 
displays significant short range order. These features are 
clear from Fig. 0] where S(k, t) is shown for e = 0.02. For 
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FIG. 6: Particle configuration snapshots at the onset of disor- 
der in the melting-freezing regime for e = 0.05, — 0.8167. 



higher values of e, the liquid- like state is more favored (as 
can be seen from the short stretches of crystalline order) 
for e = 0.06 shown in Fig. [SJ 

As is clear from Figs. El El El the time scales over 
which the crystalline order or liquid-like disorder sets 
in is considerably shorter than the persistence time of 
each of these phases. The process of ordering and disor- 
dering is better monitored through a sequence of snap- 
shots of configurations. A typical set of snapshots are 
shown in Fig. |B] The corresponding Voronoi construc- 
tion for e = 0.05 is shown in Fig. [3 ( See also Figure 5 in 
Ref.[lj|.) The extent of order in the x and y directions 
is studied using the correlation function. The nature of 
g(x) and g(y) is shown in Fig. |S]and [5] at four different 
times. Note that the extent of order in the liquid-like 
state in the direction of the drive x is significantly less 
than that in the y direction. Further the primary order- 
ing wavevectors are along i and 60° to x. Thus nearest 
neighbor distance in the y direction is y3 times larger. 
A partial understanding of why order starts to set in 
again after melting has occurred is to note that relative 
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FIG. 7: Voronoi construction from particle configurations 
showing the onset of order and of disorder for species A at 
e = 0.05, F d * = 0.8167. 




FIG. 8: The distribution function along x direction for species 
A for e = 0.05, F d * = 0.8167. 



motion disrupts primarily those structures with order- 
ing wavevector along the drift direction x, leaving some 
residual order along y as seen in Figs. and |H1 So each 
species still provides a weak periodic potential along y 
for the other species. This can induce order along x as 
well resulting in a 2d ordered state by a mechanism sim- 
ilar to the "laser induced freezing" of a 2d suspension 
of strongly interacting colloidal part icles subject to a Id 
periodic modulation [27], 0, 12H. l30| . This state persists 
for a long time, before disorder once again sets in. And it 
is in this driving force regime that one sees the melting- 
freezing cycles. Finally, we find that the two species do 
not necessarily order or disorder simultaneously. Though 
we found no clear trend in this regard, mostly, one species 
begins to order while the other is disordered. Also, for 
both species, the 'Bragg peak heights' rise fast and de- 
cline slowly, i.e., ordering takes place much faster than 
the progression of disorder. 

The structure factor can be used to obtain the dynam- 
ical phase diagram in the Fi — e plane. This is shown 
in Fig. 1101 The points shown here actually represent the 
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FIG. 9: The distribution function along y direction for species 
A for e = 0.05, F d * = 0.8167. 



value of (FJ , e) for which the crystalline phase is detected 
just before entering the bistable melt-freeze region ( re- 
gion II). The regions I and III refer respectively to the 
lower and upper sliding crystalline states. We emphasize 
here that the transition to the melt-freeze regime occurs 
over a finite range of values of the parameters and it 
is not possible to pin down with precision the point at 
which the melt-freeze alternations begin. The range of 
values of the force, F3 — F2, over which we observe the 
melt-freeze cycles increases with e and hence with Vab 
as well (Fig. I10[l . This agrees well with our observation 
that the average potential barrier that a particle has to 
negotiate during its motion in the steady sliding state 
increases with e. Further, for large values of e (e > 0.05), 
the alternations persist over a very large window of driv- 
ing forces. For such e values, we have not been able to 
detect the upper threshold F3 corresponding to the reen- 
trant crystalline state. 

We shall now discuss the mclt-frcczc cycles in more de- 
tail and explain the observed features. There is a curious 
metastability associated with the cycles: for the parame- 
ters mentioned above, a disordered configuration fails to 
nucleate even over our longest simulations if the initial 
state is chosen to be a perfectly ordered lattice. Thus, 
both the melt-free cycles and ordered sliding states dis- 
play local dynamical stability. But if we disturb this ini- 
tial perfectly ordered lattice by moving a single particle 
by, say, one lattice spacing, the melt-freeze cycles resume. 
Also note that the orientation of the triangular lattice in 
the ordered state of the cycles (Fig. Efa)) is changed by 
30° with respect to the one in the steady sliding state 
(Fig. [21 • This exchange of stabilities between the fcc-like 
(Fig. |2Jl and " layer" (Fig. Eta)) structures is known from 
experiments [3lJ and simulations ■ At low relative ve- 
locity, particles in each layer have ample time to get out 
of the way of those in the other layer while retaining on 
average the fcc-like structure. As the speed is increased, 
the structures do not have sufficient time to relax and 
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overlap is reduced by going to the "layer" structure of 
Fig. inja). However, we observe that even with such a 
shift in orientation, in the ordered part of the cycle, the 
smooth relative motion of the A and B lattices is dis- 
turbed now and then by kinks; a row moving out of step 
with adjacent rows, as marked in Fig. This leads to 
the formation of kink-like undulations transverse to the 
mean drift like a wave. At some point in time as these 
undulations build up sufficiently, the system enters a dis- 
ordered state. This state persists for a long time, before 
order sets in once again. Recall that at small e, in the 
ordered part of the cycle, the magnitude of the structure 
factor is small and correspondingly we find an enhanced 
level of fluctuations (Fig. 0} compared to that at large e 
(Fig. |3 ) . This is because the ordered states at small e 
can support a larger number of defects without making 
a transition to the disordered state, whereas as e is in- 
creased such states can not be sustained for long time. 
In fact, with increasing e, the probability of the system 
being in the ordered state crucially depends on the defect 
density, i.e., the ordered state is long-lived only when the 
number of defects is small. This can be understood by 
considering the potential landscape seen by each species. 
In the ordered state, particles of each species are nested 
in the three-fold minima formed by its nearest neighbors 
of the other species. For small e, both in the lower and 
upper part of the sliding states, the potential depth is 
shallow and creation of defects in shallow potentials does 
not cost much energy. For the same reason, the num- 
ber of such defects that the state can sustain can also 
be large which in turn implies that the extent of order 
in the crystalline part of the cycles is not significant. A 
shallow potential also allows for the ease of annealing of 
the defects as this can be accomplished by removing one 
particle from a 7-fold co-ordinated or adding one particle 
to a 5-fold coordinated site (Fig. 01 . Thus, low e situation 
allows for the ease of creation and annealing out of these 
defects as time progresses. This dynamical balance be- 
tween the creations and annihilation of the defects can be 
sustained for long stretches of time but with smaller ex- 
tent of order with concomitantly large fluctuations. This 
is precisely what is seen in Fig. 0] (Note that this pic- 
ture is also consistent with the observed feature that the 
liquid-like state at low e has a fairly high level of order 
compared to that at high e.) In contrast, with increase 
in e, the well depth increases significantly which implies 
that the crystalline order would be high as can be seen in 
Fig. 13 Moreover, the formation of defects is less favored, 
but, once formed, it cannot be easily annealed. Further 
each defect gives rise to large local restoring forces and 
the crystalline order will be terminated even when their 
number is small. 

3. For F > F3, both A and B components are once 
again well-ordered, drifting crystals. In fact, the re- 
appearance of a smooth sliding state is akin to the the 
re-entrant ordered state seen in [21, 331. 
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In fact the sliding crystalline states that we observe at 
low and high drives can be visualized as ordered states 



FIG. 10: The phase diagram of the system (100 particles of 
each species) in the e-Fa* plane, where Fd* = Flt/VAB- The 
system undergoes the stochastic 'melt-freeze' alternations in 
the region II, and is a macroscopically ordered crystal in the 
the upper and lower regions I and III respectively. 



with lanes |2Cj of single-particle width. Ref.[20j studies 
a model very similar to ours, with e = 1. In that case, 
the equilibrium state is a crystal, randomly occupied by 
each species and the driven state shows the interesting 
phenomenon of lane formation. Because e is large, the 
system tries to minimize the extent of AB interface, hence 
one gets a few broad lanes with many columns of particles 
of the same species. Since e is very small for us, we get 
many lanes of unit width. 

The results stated above are for 100 particles of each 
species. We have carried out a systematic study of this 
phenomena for 144, 169 and 256 particles of each type 
for e = 0.02. We find the same qualitative behavior as 
that for 100 particles including the range of Fd values 
for the three phases. However, for smaller systems (with 
64 particles of each species), we have not observed any 
significant decay in order for e < 0.1 (possibly due to the 
fact that the correlation length is of the order of half the 
system size in this case). 

We now discuss some natural timescales which will be 
useful in our final explanation of this phenomenon. When 
the two arrays of particles are driven through each other, 
there is a competition between two timescales, t\ the time 
to traverse one lattice spacing and t 2 the timescale of 
relaxation of a particle in the local potential well provided 
by its neighbors. For n >> r 2 , each species has ample 
time to relax to local equilibrium, while for t\ « t 2 , 
each species averages out the undulating landscape. For 
both T1/T2 » 1 and << 1, we expect and find smooth, 
orderly sliding. We therefore expect maximal effects of 
interspecies interaction where t± « t 2 (for example, for 
e = 0.05, n ~ 0.00022 and t 2 ~ 0.00019 in the 'melt- 
freeze' regime 0|) , which is what we find. This suggests 
that a detailed explanation lies in mechanisms involving 
competing timescales to which we now turn. 

The stochastic alternations of the system between 
the crystalline and liquid-like states is strongly remi- 
niscent of the phenomenon of stochastic resonance (SR) 
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22, 23, 24, 25]. To see the similarities and the differences, 
consider a prototypical example of stochastic resonance 
of a Brownian particle in a bistable potential subjected to 
a weak periodic forcing term. When half the periodicity 
of the driving force is comparable to the mean first pas- 
sage time associated with the barrier crossing, the state 
of the system switches between the two minima of the 
potential in a surprisingly regular way. Indeed, the time 
series of the position of the Brownian particle entrained 
in the two minima is very similar to S(k, t) of the particle 
model. What this potential ( or an 'effective free-energy' 
in this case) is and what modulations are, is not clear in 
the present model and needs further investigation. One 
can anticipate that the role of the bistable potential in 
SR is played by the 'effective free energy' as the system is 
a many particle system switching between the crystalline 
and liquid-like states. However, identifying a periodic 
forcing in the present context is more difficult. Thus, it 
would be useful to construct a reduced model which dis- 
plays the dominant features of the particle model. One 
important characteristic feature of SR is that the signal 
to noise ratio exhibits a maximum at an optimum value of 
the noise intensity. This aspect cannot be easily checked 
in the particle model as it involves generating very long 
time series (which would involve prohibitively large scale 
computing ). However, we recall that altering e in the 
particle model has an effect that controls the ratio of the 
residence times of the system in the ordered and disor- 
dered states. This is similar to altering the height of the 
bistable potential which in turn controls the residence 
time in the example considered. This identification fur- 
ther supports our view that the melt-freeze cycles are 
in fact stochastic resonance. We shall make this more 
concrete by introducing a reduced model which captures 
most features of the particle model. 



III. THE REDUCED MODEL 

The effects of external nonequilibrium driving condi- 
tions in an underlying first order phase transition have 
often been studied successfully by modifying, say, the 
time dependent Ginzburg-Landau equation for the dy- 
namics of the order parameter j34[. The results of the 
previous section were obtained from direct simulation of 
particle motion. In this section, we propose an under- 
standing of these results through dynamical equations 
for the appropriate order parameter fields evolving under 
the combined effect of shear and a coarse-grained free en- 
ergy. The nature of the free energy functional is usually 
determined based on the knowledge of the allowed states 
of order/disorder and a few general symmetry consider- 
ations. Here, we follow the model proposed earlier |2l| 
for studying sheared colloidal crystals [13, |3!|. Recall 
that in our simulations, at low drives, we find a smooth 
sliding crystal wherein the A lattice slides past that of 
B in a coherent fashion, and at intermediate drive val- 
ues, we observe the melt-freeze cycles. To mimic this, we 



choose an order parameter denoted by p (the Bragg peak 
intensity) , which takes on a finite value corresponding to 
the crystalline order and zero value corresponding to the 
liquid-like order. A simple form of the free energy which 
ensures the crystalline (p ^ 0) and melt phases (p = 0) 
is the Landau polynomial for a first order transition 



V(p) 



hp 3 



(5) 



The distortions produced due to the drive in the particle 
model can be represented by another variable represent- 
ing the strain (in our model, the relative phase of the 
density wave in the two sliding layers) denoted by 9. In 
the crystalline state, as distortions are small and homo- 
geneous at low drives, we take 9 to be zero for this state. 
As homogeneous distortions would mean that the suc- 
cessive minima of the crystal are equivalent, we consider 
9 to be a periodic variable with 9 = to be equivalent 
to 6 = 1. Further recall that our simulations show that 
at intermediate drives, deformation becomes inhomoge- 
neous forcing the crystalline state to melt ( although in a 
dynamical way). Again, a simple form of the free energy 
in 9 should incorporate elasticity at small 9 and yield- 
ing at large 9, ie., at small strains, V(9) is assumed to 
be quadratic in 9 and a softening term at larger strains 
( the cubic term). The coefficient of V(9) must vanish 
with p, say as p 2 , as the free energy cost of deformations 
must reduce to zero when the system is in the liquid state 
(i.e., for p = 0). The simplest general form of the free 
energy F(p, 9) of a distorted solid respecting the above 
conditions is J2j of the following form 



F{p,9) = V{p)+ l -ap 2 V{9), 
where V{9) has a similar form as V(p) 



(G) 



(7) 



The Langevin-TDGL equations for p and 9 that describe 
the dynamics of this system are 



1 dF(p,0) 
P= ~ — ^p> 



r p dp 

1 8F(p, 



T B 89 



+ n + r] B . 



(8) 
(9) 



The idea is that in the absence of any restoring forces 
for 9, 9 would be equal to il. In general, then, f2 rep- 
resents the effects of relative sliding of the layers, at a 
rate determined by competition between and ^jf-. If 
is too small, 9 will get stuck at a finite value in the 
absence of noise. r g s (q = p,9) are the kinetic co- 
efficients, 77 g 's represent Gaussian delta-correlated noise 
components whose variances are related to T q and tem- 
perature through the fluctuation-dissipation relation. We 
ignore possible additional non-equilibrium noise sources. 
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The equations for our system are those for an over- 
damped particle in an nonsymmetric double well poten- 
tial F(p, 9), driven along the angular coordinate [2a]. F° r 
zero strain (9 = 0), the system relaxes in either of the 
two minima corresponding to the liquid minimum 



Pi =0,9 = 0, 
or the crystalline minimum 



Pc 



2ci 




4aiCi 

7 2 

01 



0, 



(10) 



(11) 



depending on which is the locally stable state. 

We choose the parameter values (oi, b\ and c±) such 
that the crystalline minimum p c is the more favorable 
state at zero drive and the potential barrier between the 
two minima, V(pi) — V{p c ), is appropriate. In the pres- 
ence of noise ( whose strength can be appropriately cho- 
sen) the system equilibrates with the respective popula- 
tions determined by noise strength and the relative well 
depths. If 9 and p are macroscopic (infinite system-size) 
averages, there should be no noise in the equations. In 
practice, presumably, shear melting and the cycles take 
place over a finite correlated domain (whose size we do 
not know). We are therefore justified in using a noisy 
evolution equation. 

When the drive f2 is switched on, this scenario is al- 
tered and the populations of each of these wells now 
evolve in time depending upon the competing time scales 
of relaxation and applied shear rate. In this case, it is 
better to consider the fixed points of the noise free case 
of Eq. [S] and ED The two attractive fixed points to which 
the system relaxes can now be identified with crystalline 
and liquid-like order. The repulsive fixed point deter- 
mines a saddle type of maximum. These will depend on 
the shear rate fi. The barrier height between the stable 
fixed points and unstable fixed point determine the barri- 
ers that the system has to surmount. These now depend 
on 0. We find that the value of the free energy at the 
liquid like minima and the saddle do not change signifi- 
cantly, only the crystalline (distorted) minimum changes 
as a function of Q given by 



2ci 




4(ai + aV{9 min ))c 1 
oi 



(12) 



Note that Eq. ED determines a critical value of 
f2 = f2 c = 0.5ap 2 ^-\ max . For values of < f2 c , the 
barrier is high. In such a situation, in the presence of 
noise, the transitions would be rare. But, on increas- 
ing f2 beyond the critical value the "free energy" of 
crystalline state becomes comparable with that of the 
liquid minimum. Under these conditions, noise assisted 
transitions to the liquid state occur. More importantly, 
in this regime, as 9 is itself changing as a function of 
time, the minima in (p, 9) is slowly modulated under 
action of the drive f2 and the relative stability changes 



as a function of time. When the time scale imposed by 
Q. is small enough to allow the system to make inter-well 
transitions and when the time scale of the induced 
periodicity is a ppr oximately equal to the Kramers 
escape time |37l |38( under the influence of noise, one 
expects transitions between the crystalline and liquid 
like states in a range of values of beyond fl c . As 
a result the system can undergo stochastic transitions 
between the two metastable states which in turn can 
lead to comparable lifetimes. 

As this is a driven system, a reasonable criterion 
for studying the occupancy of the system is to cal- 
culate the marginal probability distribution function 
P{p) = I P{p,9)d9 (i.e., the probability of the order 
parameter having the value p, independent of the value 
of the strain field 9). In the following section we show 
that the external drive causes the system to sample both 
minima or stay mainly in one of the minima depending 
on the value of the drive Q, noise strength, coupling 
constant a starting from an initial crystalline state. 



A. Results of the Reduced Model 

We study the time series and the probability distribu- 
tion of the system by discretizing the Langevin equations 
in time. The integration scheme used is the fourth order 
Runge-Kutta with a fixed time step of 0.001. After dis- 
carding transients (~ 10 5 time steps) the time evolution 
of the system for the next 8 x 10 9 time steps is moni- 
tored. In figures ITT1 and IT31 we have shown a timeseries 
stretch from 6 x 10 9 to 6.5 x 10 9 timesteps. Noise cor- 
responding to p and 9 are drawn from Gaussian white 
noise distributions with zero mean and unit variance. 
For studying the time evolution of the system, we have 
chosen a noise strength D = 7 x 10~ 4 for both rj p and 
i]0. ( We shall study the influence of noise on the signal 
to noise ratio later.) For the numerical work reported 
here, the parameters values in the free energy used are 
ai = 0.85, 6i = 5.8, c\ = 8.0, and for the strain field are 
a 2 = 1.3644, b 2 = 8.7105 and c 2 = 13.6740. We study 
the time evolution of the system as a function of the drive 
fl (for a fixed coupling parameter a) and as a function of 
a for a fixed drive. We observe the following behavior: 

a. As a function of the drive fl : 

Keeping a fixed at an optimum value (a = 0.17), the 
drive has two thresholds: 

i. For £1 < fii and £1 > O2, we find the system always 
resides in the crystalline minimum (p 7^ 0). The time se- 
ries of pit) is shown in Fig.lllb.c for a typical set of values 
fl = 0.03 and fl = 0.3 respectively. The corresponding 
probability distribution is peaked around the crystalline 
minimum, as can be seen from Fig. 112b . c. 

ii. In an intermediate window of driving forces < 
£1 < Q 2 , the system stochastically alternates between the 
liquid (p = 0) and the crystalline minimum. The time 
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FIG. 11: The time series of the order parameter in the low 
(top panel), high (bottom panel) and intermediate (middle 
panel) driving force regimes. We have chosen a noise strength 
D = 7 x 1CP 4 and coupling parameter a = 0.17. 
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FIG. 13: The time series of the order parameter in the low 
(top panel), high (bottom panel) and intermediate (middle 
panel) values of the coupling parameter a for noise strength 
D = 7 x 10~ 4 and driving force Q = 0.08. 
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FIG. 12: The marginal probability distribution of the order 
parameter in the low (top panel), high (bottom panel) and 
intermediate (middle panel) driving force regimes for the same 
parameter values as in Fig.ll. 



series corresponding to a typical value of fl — 0.08 is 
shown in Fig. Illb . The probability distribution has two 
peaks corresponding to these two minima. The time se- 
ries of the order parameter (Fig. ITTb ) also shows that the 
persistence time of these two states are comparable (for 
optimum values of a) very much like the time series of 
a system undergoing stochastic resonance. These results 
are similar to the results of the particle model where F£ 
was varied for a fixed value of e. 

b. As a function of the coupling a : 

Here, we have kept the drive fl at intermediate values. 



We find that the crystalline minimum is favored over the 
liquid one for low values of a, whereas for high values 
the liquid minimum dominates. At intermediate values 
of a, the system spends comparable durations in each 
state. As we increase a from small values, one finds that 
the system tends to spend increasingly more time in the 
liquid state and, eventually, at large values of a the liq- 
uid state is the preferred state. This is shown in Fig. IT3l 
for three typical values of a. (The numerical results are 
for £1 = 0.08 for which we observe the most prominent 
stochastic switching of the order parameter values be- 
tween p ^ and p = 0.) This is also reflected in the 
probability distribution which has a more pronounced 
peak at the crystalline minimum at small a and at the 
liquid minimum at large a. For intermediate values of a 
we find a bimodal distribution. Figure lHI shows the prob- 
ability distributions for the parameter values of Fig. 1131 
This behavior is similar to the results obtained in the 
particle model by varying e keeping F£ fixed. 

As mentioned, one dominant feature of SR is the en- 
hancement of the signal to noise ratio (SNR) at an op- 
timum value of the noise intensity. In order to check 
this, we have carried out long runs of the order of 10 10 
time steps. The power spectral density (PSD) of the 
time series has been calculated for various values of the 
noise intensity D — 6 x 10~ 4 to 1.2 x 10~ 3 . It exhibits 
strong peaks at all integral values of the fundamental 
(unlike the symmetric bistable potential where only the 
odd harmonics are seen) due to the absence of any sym- 
metry in F{p,ff). A plot of this is shown in Fig. ^| 
for a typical value of D = 7.0 x 10~ 4 . The signal to 
noise ratio calculated from the power spectrum using the 
first peak for various values of D is shown in Fig. 1161 
(Here we have used the conventional SNR defined by 
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FIG. 14: The marginal probability distribution of the order 
parameter in the low (top panel), high (bottom panel) and 
intermediate (middle panel) a, for the same parameter values 
as in Fig. 13. 
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FIG. 15: The power spectrum of the time series of the order 
parameter p for a = 0.17, Q = 0.08 and D — 7.0 x 10" 4 . 



SNR = lt)log[S signa i(uj) / S noise {oj)].) As is clear from 
the figure, the maximum enhancement of the SNR is 
found to be around D ~ 7.0 x 10 -4 . 

As discussed above, we note that fi and a of the re- 
duced model take the roles of the drive and the in- 
terspecies interaction strength e in the particle model re- 
spectively. To make the parallel between these models 
more concrete, we have constructed the dynamical phase 
diagram in the a — Q plane shown in Fig. 1171 ( In con- 
structing this diagram, we have taken the system to be 
a liquid state if it spends less than 2% of the time in the 
crystalline minimum, and correspondingly for the crys- 
tal.) Region I refers to the crystalline phase and region 
III, the re-entrant crystalline phase. The melt-freeze cy- 
cles where the system alternates between crystal and liq- 
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FIG. 16: The signal to noise ratio for the reduced model for 
parameter values a — 0.17, SI = 0.08. 
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FIG. 17: The phase diagram of the reduced model in the a-Q, 
plane. The system is crystalline in region I and III represents 
the re-entrant solid. In region // we observe the 'melt-freeze 
cycles'. For high values, liquid-like region IV is seen again. 



uid is shown as region II. In region II, for low values of 
the coupling parameter a, the persistence of the ordered 
state is more than that for the disordered state, and de- 
creases as a is increased, eventually giving way to the 
liquid- like region IV for moderate values of fl. It is clear 
that this diagram is remarkably similar to the phase dia- 
gram of the particle model shown in Fig. 1101 except for 
the region IV which could not be detected in the particle 
model within the longest runs carried out. 

The physical picture of the reduced model is clear. As 
the initial state at zero shear rate is a crystalline state, 
by continuity arguments one should expect that at low 
shear rates the system should find the crystalline mini- 
mum favorable. At intermediate range of shear rates, the 
system develops another minimum at zero value of the or- 
der parameter p corresponding to the melt state making 
the system bistable. When the shear rate is close to ( 
and larger than) the critical value f2 c , under the action 
of noise, the systems makes transitions from the crys- 
talline minimum to the liquid minimum and vice versa. 
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FIG. 18: AF C and AF t as a function of Q. 



However, since the strain 9 evolves in time the system 
experiences an additional periodic modulation. We note 
here, that this periodicity is not equal to l/f2 as the the 
strain variable moves on the V(9) surface. Typical values 
of the induced periodicity estimated from deterministic 
version of Eqs. IHlandElfor the optimum range of drive 
values is of the order of 10 5 time units. When the induced 
periodicity is twice the Kramers escape rate [32j , the sys- 
tem alternates between the two minima. Further, we note 
that as the two wells are not symmetric, in the presence 
of noise, the mean first passage times associated with the 
two wells will be different. Indeed, we find that as a func- 
tion of O, there is a range of f2 values (0.065 < £1 < 0.12) 
where the barrier between the crystalline minimum and 
maximum of the free energy F sadd i e - F crysta i = AF C , 
is much smaller than that between the liquid like mini- 
mum and the maximum, F sadcUe - F liquid = AF ; . A plot 
of AF C and AF; is shown in Fig. EH It is clear that 
AFi is nearly constant as a function of fl while AF C goes 
through a minimum in the range of fJ = 0.065 to 0.12. A 
simple order of magnitude calculation gives the Kramers 
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2tt 



-exp - (AF/D) ~ 10" 



which matches with the frequency range of the first har- 
monic seen in the PSD (see Fig. 115(1. We note here that 
the first passage time for the system to cross the bar- 
rier from the liquid state is much too large as the barrier 
AFi ~ 0.01. Thus, it is clear that the reduced model 
reproduces, qualitatively, most features of the particle 
model. The time series of S(k, t) in the particle model is 
very similar to that of p(t) which we have demonstrated 
has all the features of stochastic resonance. The regime 
showing stochastic resonance is more pronounced for op- 
timum values of f2 and a. In fact here, as in our particle 
model simulations, in the melt-freeze regime, for low val- 
ues of the coupling parameter, the crystalline state is fa- 
vored and for high values of the coupling parameter, the 
liquid state is favored, while in an intermediate regime 
they enter roughly equally. The similarity of the reduced 
model with the particle model is well summarized by the 



phase diagram of the model in the — a plane which is 
similar to that of the particle model in the F*[ — e plane. 



IV. DISCUSSION 

In summary, we have studied the nonequilibrium sta- 
tistical behavior of two adjacent monolayers sheared past 
each other, using Brownian dynamics simulations. For 
low and high driving forces, we obtain macroscopically 
ordered, steadily drifting states. In a suitable range of 
driving rates we see that the system switches between 
crystalline and liquid like states. The residence times 
are nearly equal in a intermediate range of values of the 
inter-layer coupling. As we have seen, the interlayer cou- 
pling essentially determines the barrier each particle has 
to surmount in order keep pace with the applied drive. 
Our simulations show that the switching between the 
crystalline and liquid like states maintains the spatial 
coherence ( or the lack of it in the liquid state) and thus 
is an example of cooperative stochastic resonance. Al- 
though there is no external imposed time-periodic po- 
tential present in our system, its role is played by the 
drive which shears the two layers past each other. Thus, 
each moving layer provides a time-varying potential. To 
make the connection to SR more concrete, we have in- 
troduced a reduced model for the dynamics of amplitude 
and phase of the crystalline order parameter |2l| which 
mimics all features of the particle model apart from show- 
ing the main feature of SR namely the enhancement of 
signal to noise ratio for optimum values of noise strength. 
This relation between the particle model and the reduced 
model makes a convincing case that the melt-freeze cy- 
cles observed in the former are indeed a manifestation of 
stochastic resonance of a spatially extended noisy inter- 
acting system subjected to a constant drive. 

A few comments may be in order on the dynamics of 
the particle model. Although our model falls broadly 
into the class of shear driven systems which are known 
to display stick-slip behavior 0, fl3l l39| in most such 
cases the system spends considerable time in the stuck 
phase and very little time in the slip phase. Experiments 
on confined fluids of few monolayers (a situation superfi- 
cially similar to our model) and the corresponding molec- 
ular dynamics simulations [l2l IT^ | show that the system 
also displays melt-freeze cycles with the crystalline phase 
persisting significantly much longer than the melt phase. 
From this point of view, the persistence dynamics we ob- 
serve with nearly equal residence times in the crystalline 
state and the liquid-like states comes as a surprise and 
is clearly not conventional stick-slip. The reduced model 
has helped us to elucidate the connection to stochastic 
resonance apart from the similarity of the phase diagrams 
of the two models. However, it is clear that not all as- 
pects of the particle model are mimicked by the reduced 
model. For instance, the feature observed for small in- 
terlayer coupling e, namely the smaller extent of order as 
seen by the small values of S(k, t) in the crystalline state 
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( compared to larger values of e) and larger fluctuations is 
not seen in the reduced model. As explained earlier, this 
feature is a many body effect appearing in a low barrier 
situation and therefore cannot be explained on the basis 
of the SR features of the reduced model valid only in the 
Kramers high barrier limit (even if one were to include 
appropriate spatial degrees of freedom). 

Finally, let us comment on the experiments which can 
verify our simulations. We expect that this phenomenon 
should arise in adjacent crystal planes of sheared col- 
loidal crystals. To study this effect in bulk sheared col- 
loidal crystals would require not conventional scattering 
probes but something that focuses on an adjacent pair 
of crystal planes . Alternately we could look at two 
solid surfaces patterned with ordered copolymer 0, 0] 
or colloidal monolayers under relative shear. The colloids 



or copolymer adsorbed to the two surfaces must be of two 
different kinds, each having more affinity towards one of 
the surfaces. Another possibility would be to study elec- 
trokinctic motion of charged 2d confined binary colloids 
in the presence of a constant external electric field, gen- 
eralizing the ideas of (2£j • Centre of mass measurements 
would be the same as that made for two species being 
driven in opposite directions. 
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